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Dynamics of Domain Growth in Self- Assembled Fluid Vesicles 
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The dynamics of phase separation in multi-component bilayer fluid vesicles is investigated by 
means of large-scale dissipative particle dynamics. The model explicitly accounts for solvent par- 
ticles, thereby allowing for the very first numerical investigation of the effects of hydrodynamics 
and area-to-volume constraints. We observed regimes corresponding to coalescence of flat patches, 
budding and vesiculation, and coalescence of caps. We point out that the area-to-volume constraint 
has a strong influence on crossovers between these regimes. 
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Apart from partitioning the inner and outer environ- 
ments of the cell, biomembranes also act as support 
for complex and specialized molecular machinery, crucial 
for variousphysiological functions and trans-membrane 
transport U| . It is believed that membranes maintain in- 
plane compositional organization that is essential for its 
function. Though biological membranes are very complex 
and can in general be far from equilibrium, knowledge of 
the equilibrium properties of simple model membranes 
will be extremely useful in providing understanding of 
this molecular machinery. Such a study should also be 
essential for developing new applications involving lipo- 
somes. This has recently resulted in a surge of theoretical 
and experimental investigations on domain formation in 
multicomponent vesicles jllllllllllllllllIS- 

The dynamics of in-plane demixing of multicomponent 
vesicles into coexisting phases is richer than its counter- 
part in bulk systems. There are several reasons for this: 
(i) the phase separation process is strongly coupled to 
the shape dynamics of the vesicle, (ii) the viscosities of 
the lipid bilayer and that of the embedding solvent are 
different, and (iii) the membrane can be impermeable to 
solvent, resulting in a constraint on the vesicle area-to- 
volume ratio. The purpose of this letter is to investigate 
the effect of the above on the dynamics of multicompo- 
nent fluid vesicles. 

Phase separation of vesicles following a quench to 
the two-phase region has previously been considered by 
means of a generalized time-dependent Ginzburg-Landau 
model H, and dynamic triangulation Monte Carlo 
model 0, la Q- The later study showed a marked de- 
parture of the phase separation in multicomponent fluid 
vesicles from their Euclidean counterparts as a result of 
coupling between curvature and composition. In particu- 
lar, it was shown that at intermediate times, the familiar 
labyrinth-like spinodal network breaks up into isolated 
domains 5, 6]. In the case of a tensionless closed vesicle, 
at late time these domains reshape into buds connected 
to the parent vesicle by very narrow necks Q- Further 
domain growth proceeds via Brownian motion of these 
buds and their coalescence. 



Recent advances in experimental techniques, like the tide i: (i) a conservative force, ¥>■ , (ii) a dissipative 



two photon fluorescence and confocal microscopy, has 
made it possible to study phase separation on fluid vesi- 
cles l^, S ll^- These experiments reported structures 
with many domains which arc more akin to caps than 
fully developed buds. A natural question then is: Does 
coarsening in multicomponent fluid vesicles proceed via 
a kinetic pathway similar to that predicted by the recent 
simulations? 

Previous studies of the dynamics of phase separation 
in fluid membranes did not take into account the fol- 
lowing important features of a real lipid vesicle: (i) the 
presence of solvent, (ii) the constraint of area-to- volume 
ratio and (ii) the freedom to vesiculate. We report here 
the first large scale simulation of phase separation dy- 
namics of a fluid vesicle model that accounts for these 
features. This study, based on dissipative particle dy- 
namics (DPD), finds that the dynamics at all times is 
affected by the presence of the solvent, and that the late 
time shape of the vesicle corresponds to that of a surface 
decorated with caps. 

Within the DPD approach [ij, the vesicle is formed 
from the self-assembly of individual particles in an ex- 
plicit solvent. The model parameters are chosen such 
that the membrane is impermeable to the solvent, thus 
allowing us to investigate the effect of conservation of 
both vesicle's area and encapsulated volume. A lipid 
particle is modeled as a simple flexible amphiphilic chain 
built with four DPD particles; one hydrophilic particle, 
simulating the lipid head group (h) , and three hydropho- 
bic particles, simulating the lipid hydrophobic part (t). 
The heads of the two types of lipids are denoted by Ha 
or IiB and their tails are denoted by Ia or Ib- Water 
particles are denoted by W . We focus on the case where 
interactions are symmetric under exchange of A and B 
components, thus ensuring no explicit coupling between 
local composition and local curvature. 

The position and velocity of each particle are denoted 
by (r^ , Vi ) . Their time evolution is governed by New- 
ton's equations of motion [l3|. There are three types of 
pairwise additive forces acting on a particle j by a par- 
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FIG. 1: Snapshots of a phase separating vesicle in case-I at 
(a) i = 100 r and (b) 2000 r. (c) is a slice taken at 5000 r. 
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force, Fj- and (c) a random force, F^, . The conser- 
vative force between two particles i and j is given by 
py — —aijUj{rij)nij where we choose uj{r) = 1 — r/vc 
for r < Tc and Lo{r) =0 for r > r^ such that the forces 
are soft and repulsive. Here r.y — Yj — r^, riy — r,y/|ry | 
and Tfj is the cutoff radius for the interaction. The hy- 
drophobic and hydrophilic interactions emerge from the 
relative interaction strengths a^. With this, the param- 
eters chosen for the simulation are 
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where e sets the energy scale. The interaction parameters 
are chosen such that the amphiphiles self-assemble into 
bilayers and the two types of lipids are in the strong seg- 
regation regime. A lipid particle integrity is ensured via 
a harmonic interaction between consecutive monomers 
within a chain given by F^ -:;_;^ = -C(l - ri^i+i/b)ni^i+i, 
where C is a positive spring constant, and 6 is a preferred 
bond length. We use b = 0.45rc and C = lOOe. The dissi- 
pative force originating from "friction" between particles 
is given by, Flf' = -7^ w2(ry)(ny • Vy)ny where 7^ is a 
friction parameter and Vy ~ Vj — v^ . The random force 



is given by F, 
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-aijU!{rij)QjiAt) ^/^n^, where. At 



is the time step in units of t = (mr^/e)^' ^ with m being 
the mass of a single DPD particle. All particles in our 
simulation have the same mass. Here, dj is a symmetric 
random variable with zero mean and unit variance, un- 
correlated for different pairs of particles and at different 
time. The dissipative and random forces are related to 
each other through the fluctuation dissipation theorem, 
leading to the relation jij — afj/2kQT. Unless otherwise 
specified, we used a fixed value cry — a for all interacting 
pairs. The pairwise nature of the dissipative and random 
forces ensures local momentum conservation leading to 
correct long-range hydrodynamics 13]. 

In our simulation the following parameters were used 
a = 3.0(e3TO/r2)i/4^ At = 0.05t and /csT = e. With 
the parameters chosen, the estimated bending modu- 
lus, K ^ lO/c^T, and line tension, A ~ 10~^^J/^m, are 
in reasonable agreement with estimates for lipid mem- 
branes [ij|. We use 16,000 lipid chains in a simulation 



cle, when it is not deflated, is about 138,400. The total 
number of DPD particles corresponds to 1,536,000 ^\w\- 
Bilayers and vesicles can be self-assembled in a DPD 
simulation [lfl[i3- In order to save computer time, we 
prepare our vesicles starting from an almost closed con- 
figuration, composed of a single lipid component, and 
allow it to equilibrate. This approach, allows for the 
equilibration of the lipid surface density within the two 
monolayers through the edge of the membrane. We find 
that the vesicle closes within about 50,000 DPD steps. 
Within this model parameters, we found that the num- 
ber of solvent particles inside the vesicle and the num- 
ber of lipid particles in the inner and outer leaves of the 
vesicle remain constant. This implies that the vesicle is 
impermeable to the solvent and flip-fiop events are rare. 
The vesicle prepared in this manner is found to be nearly 
spherical with about 140, 000 solvent particles inside it. 
We refer to configurations with this inner volume as case- 

I. A vesicle with excess area can then be prepared by 
transferring a certain amount of solvent particles from 
its inner to outer regions. We refer to as case-H the equi- 
librated configurations obtained after 20% of the inner 
solvent particles are transfered. We found that in case- 

II, as well as in case-I, the solvent density is the same 
inside and outside the vesicle. 

To study the phase separation, we randomly label each 
lipid particle in the bilayer as type A or i3, such that the 
relative composition in each layer is the same. Coars- 
ening is subsequently monitored by computing the dis- 
tribution of cluster sizes of the minority component and 
the total interface length between the domains. We will 
focus on the case of off-critical quenches with the volume 
fraction of the B-component equal to 0.3. 

Snapshots of the time evolution of two component vesi- 
cles in cases I and II are shown in Figs. ^ and |21 respec- 
tively. As shown in Fig. ^ (c), domains on either side 
of the membrane are in register. This inter-layer align- 
ment since early stages of the phase separation dynamics. 
The net interface length is shown in Fig. O During rela- 
tively early times, in both cases I and II, domains are fiat 
circular patches. The net interface length, L(t), during 
this regime is independent of the area-to-volume ratio, 
and has the form L(t) ^ t^", with a ~ 0.3, as shown 
in Fig. 121 In this regime, the number of clusters on the 
vesicle, N^.^ is also found to obey a power law, Nc ^ t~^ , 
with P « 2/3. This is depicted in Fig.H 

Since the amount of A and B components on the mem- 
brane are conserved, we have Nc{t) '^ AB/R'^{t), where 
R is the average domain size and Ab is the total area 
occupied by the B patches. Since the patches are fiat 
R{t) ^ AB/L{t). Thus one expects f3 = 2a. The data 
shown in Figs. Oand^lare consistent with this relation. If 
the diffusion coefficient of disks moving on the membrane 
depends on the domain size as Dd ^ l/R{t) then coales- 
cence of domains leads to a = 1/3 and /3 = 2/3 [l8| . 
These exponents can also be the result of evaporation- 







FIG. 2: Snapshots of a phase separating vesicle in case-II. Snapshots from left to right are at t 
4000 r, respectively. Analysis show that the detached objects are full vesicles. 
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FIG. 3: Net interface length as a function of time in case-I 
(blue curve) and case II (red curve). The green curve corre- 
sponds to case II with low line tension (see text for details of 
corresponding interactions). The dashed lines and the dot- 
dashed line have slopes of 0.3 and 4/9, respectively. In the 
inset, the number of vesiculated domains is shown as a func- 
tion of time for case-II with high line tension (corresponding 
to the red curve). The arrows point to the time regime during 
which budding and vesiculation occur in case II with high line 
tension. 



condensation mechanism, as predicted by Lifshitz and 
Slyozov [13 ■ However, we could monitor many events 
of patch coalescence in our simulation. We thus believe 
that the mechanism of domain growth during this regime 
is coalescence of patches ^y ■ 

After about t = 400r, the dynamics in cases I and 
II clearly departs from each other, as shown in Fig. |31 
In case-I, although there are some caps formed, coarsen- 
ing proceeds mainly through coalescence of flat circular 
patches (see Fig. QJ. Further capping is suppressed by 
the lateral tension on the membrane due to volume con- 
straint. Similar capping-induced tension was also seen 
in earlier Monte Carlo simulations @|. In case-II, where 
there is more excess area, domains reshape into caps. 
The presence of excess area allows a fraction of the caps 
to further reshape into buds which then vesiculate. The 



buds vesiculate during a short time implying that the 
energy barriers involved are very small. We confirm this 
by performing simulations on a vesicle with large excess 
area having a single B component domain occupying 12% 
of the total area. Once the bud is formed, it is found to 
vesiculate within a time period of lOr. This budding and 
vesiculation results in a marked decrease in the net in- 
terface length L{t) over a short period of time as shown 
in Fig. 13 Pinching mechanism of vesicles was recently 
studied in detail using DPD simulations \V^ . 

All detached vesicles are composed of the B compo- 
nent. Budding occurs towards both inner and outer re- 
gions of the vesicles, although most buds are found in the 
outer region. Inward budding is due to the bilayer na- 
ture of the vesicles and cannot be observed if the vesicle 
is modeled by a single surface. In the absence of explicit 
coupling between mean curvature and composition, in- 
ward vesiculation of a B-domain results when the area of 
the domain on the inner monolayer exceeds its area on 
the outer monolayer. 

Vesiculation results in a reduction of area-to-volume 
ratio of the parent vesicle leading it to acquire a more 
spherical shape. The lateral tension resulting from this 
prevents further increase in curvature of caps. Coarsen- 
ing now proceeds via coalescence of these caps. Since 
patch and cap coalescence, in a viscous medium, lead to 
the same exponent, (3, we do not see any marked change 
in the behavior of the curve shown in Fig. 2|for case-II. 
During this regime, the net interface length L ^ i°-^, 
as shown in Fig. |3| The exponent a ~ 0.4 is between 
that due to coalescence of flat patches, i.e. a — 1/3, 
and that due to coalescence of buds having an interface 
length independent of the bud size, i.e. a = 2/3 [3. As 
can be seen in Fig. 2, domains remaining on the vesi- 
cle typically have a shape closer to a hemispherical cap 
than a complete bud. In this case, the average inter- 
face length of a cap, Ic '^ (ko/A)^/"^, where a is the cap 
area. In the absence of further vesiculation, Nc a is a 
constant of time. We thus obtain for the net interface 
length L ~ NJc ~ t~^/^ , consistent with our numerical 
results shown in Fig. |3| 

Now, let us discuss the effect of line tension on do- 
main growth in case-II. For a given _B-cap of area a, sur- 
rounded by the ^-phase, the bending energy of the cap 




FIG. 4: Number of clusters on the vesicle as a function of 
time. Curve colors correspond to those in Fig. |3| The slope 
of the dashed line is 2/3. 




FIG. 5: Snapshots sequence for case-II, with line tension lower 
that in Fig|21 Snapshots from left to right correspond to t = 
1200, 3700, and 7500 r, respectively. 



scales with its radius of curvature, r, as na/r^ . Since the 
perimeter of the domain decreases with decreasing r, the 
cap adopts a shape with finite curvature. The domain 
area beyond which capping occurs is around 47rK^/A^. 
Let ii be the time required to reach this capping regime. 
Beyond ti, domain growth proceeds via cap coalescence. 
In the presence of enough excess area, as in case-II, this 



coalescence may lead to caps with interface length com- 
parable to the thickness of the bilayer. Vesiculation may 
then proceed \u\ . 

For a A, and k, there is an area of the cap around 
167rK^/A^ beyond which a complete bud is formed, pro- 
vided that there is enough excess area. Let i2 be the 
average time necessary to reach this area. During the 
period between ti and ^2, growth should be mediated by 
coalescence of caps. In case-II, discussed above, t2 — ii 
is too short to numerically detect this regime. We there- 
fore performed another set of simulations with a reduced 
value of the line tension, achieved by choosing the in- 
teraction parameter Oy for t^ts and hAhs pairs to be 
50e/rc, Corresponding snapshots are shown in Fig.|Sl As 
can be seen from this figure, vesiculation does not occur 
even at very late times. However, well defined caps with 
finite curvature are prominent. The net interface length 
shows a clear crossover between the patch coalescence 
regime and caps coalescence regime, as shown in Fig. |31 

In summary, we presented here the very first study of 
phase separation dynamics of self-assembled bilayer fluid 
vesicles with hydrodynamic effects. We found rich dy- 
namics with crossovers that depend strongly on area-to- 
volume ratio and line tension between the two coexisting 
phases. In particular, wc found that vesiculation, when it 
happens, occurs during a short period of time. Another 
feature that was not possible to account for in earlier sim- 
ulations but has important consequences is the volume 
constraint. As a result of tension induced by this con- 
straint, the main vesicle eventually acquires a spherical 
shape decorated with caps, irrespective of initial area- 
to-volume ratio, a result seemingly born out in recent 
experiments |8ll9l[l0l| 
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